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ABSTRACT 

We use an extremely large volume (2.4/i _3 Gpc 3 ), high resolution N-body simulation 
to measure the higher order clustering of dark matter haloes as a function of mass 
and internal structure. As a result of the large simulation volume and the use of a 
novel "cross-moment" counts-in-cells technique which suppresses discreteness noise, 
we are able to measure the clustering of haloes corresponding to rarer peaks than was 
possible in previous studies; the rarest haloes for which we measure the variance are 
100 times more clustered than the dark matter. We are able to extract, for the first 
time, halo bias parameters from linear up to fourth order. For all orders measured, we 
find that the bias parameters are a strong function of mass for haloes more massive 
than the characteristic mass M* . Currently, no theoretical model is able to reproduce 
this mass dependence closely. We find that the bias parameters also depend on the 
internal structure of the halo up to fourth order. For haloes more massive than M*, 
we find that the more concentrated haloes are more weakly clustered than the less 
concentrated ones. We see no dependence of clustering on concentration for haloes 
with masses M < M*; this is contrary to the trend reported in the literature when 
segregating haloes by their formation time. Our results are insensitive to whether 
haloes are labelled by the total mass returned by the friends-of-friends group finder or 
by the mass of the most massive substructure. This implies that our conclusions are 
not an artefact of the particular choice of group finding algorithm. Our results will 
provide important input to theoretical models of galaxy clustering. 
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1 INTRODUCTION 

The spatial distribution of dark matter haloes is not as sim- 
ple as was once suspected. In the standard theoretical model 
for the abundance and distribution of haloes, the cluster- 
ing strength of haloes is predicted to be a function of mass 
alone, with more massive haloes displaying stronger clus- 



terin g (e.g. lKaisedll984l ; ICole fc Kaisenll989l ; iMo fc Whitel 
1 19961 ) . However, recent numerical simulations of hierarchi- 
cal cosmologies, by covering larger volumes with ever im- 
proving mass resolution, have been able to reveal subtle 
dependences of halo clustering on other properties such as 
formation redshift, the internal structure of the halo and 



its spin dGao, SorinEel. & Whitel 


2005 


;IWechsler et al.l 


2006 


Harker et alj 20061; Bett et al.l 


2007 


; Wetzel et al.l 


2007 


Jine, Suto, & Mo 2007; EsDino-Brioncs et al. 2007h. 



The dependence of halo clustering on a second param- 
eter in addition to mass is generally referred to as assembly 
bias. However, the nature of the trend in clustering strength 
recovered depends upon the choice of property used to clas- 
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sify haloes of a given mass. Early simulation work failed 
to uncover a convincing assembly bias signal, as a result of 
insufficient volume and mass resolution, which meant that 
halo clustering could be measured f or only a narrow range 
of mass and with limited statistics dLemson fc K auffmann 
1 19991 ; iPercival et ai"1l2003l ; ISheth fc Tormenll2004l ) The first 
clear indication of a dependence of halo clustering on a 
second property was uncovered by iGao. Springel. fc White] 
|2005l ). These authors reported that low mass haloes which 
form early are more clustered than haloes of the same mass 
which form later on. No effect was seen for massive haloes. 
Wcchslcr et all l|2006h were able to confirm this result but 
also found that halo clustering depends on the density pro- 
file of th e halo, as characterized by the concentration pa- 
rameter l|Navarro. Frenk. fc Whitel Il997l 'l . The sense of the 
dependence of clustering strength on concentration changes 
with mass. Wechsler et al. found that massive haloes showed 
a dependence of clustering strength on concentration, with 
low concentratio n haloes being the more strongly clustered 
(as confirmed bylGao fc Whitefl2007l . |jing. Suto. fc Moll2007l 
and I Wetzel et al.l 120071 ). This trend of clustering strength 
with concentration is reversed for low mass haloes. Al- 
though formation time and concentration are correlated (e.g. 
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iNeto et all [20071 ), their impact on the clustering of haloes 
does not follow trivially from this correlation, suggesting 
that some other parameter may be m ore fundamental (as 
argued bv lCroton. Gao. fc Whitell2007l ). 

Previous studies of assembly bias have focused exclu- 
sively on the linear bias parameter, which relates the two- 
point correlations of haloes and dark matter. Measurements 
from local surveys have shown that galaxies have signifi- 
cant higher order correlation functions and that the spa- 
tial distribution of galaxies and haloes is not fully described 
by two-point statistics (e.g. Baugh et al. 20041 ; ICroton et al.l 
12004 iNichol et al.ll2006l ; iFrith et al.l 12006ft . With large sur- 
veys planned at higher redshifts, there is a clear need for 
accurate models of the higher order clustering of dark mat- 
ter haloes, and to establish whether or not the higher order 
bias parameters depend on other properties in addition to 
mass. 

In this paper we measure the higher order bias pa- 
rameters of dark matter haloes using a simulation which 
covers a volume more than an order of magnitude larger 
than the run analyzed by Gao and collaborators. We use a 
novel approach to estimate the higher order correlation func- 
tions of dark matter haloes. Our method builds upon the 
cross-cor relation technique advocated for two-po i nt cor rela- 
tions bvlJing Suto. fc Mol |2007MGao fc White! (|2007h and 
ISmith et al.l ( 2007 ). By considering fluctuations in the den- 
sity of haloes and dark matter within the same smoothing 
window, we can suppress discreteness noise in our measure- 
ments. This improved clustering estimator, which uses the 
counts-in-cells method, when coupled with the large volume 
of our simulation, allows us to recover the bias parameters 
from linear to fourth order, and to study the dependence of 
these parameters on the halo concentration. 

In Section 2, we give the theoretical background to the 
counts-in-cells technique we use to estimate higher order 
clustering and explain how the clustering of haloes relates to 
the underlying dark matter at different orders. We also intro- 
duce the numerical simulations in that section. We present 
our results in Section 3 and a summary and discussion in 
Section 4. 



2 THEORETICAL BACKGROUND AND 
METHOD 

In this Section we give the theoretical background to the 
measurements presented in Section 3. We estimate the clus- 
tering of haloes and dark matter using a counts-in-cells ap- 
proach. An overview of this method is given in §2.1, in which 
we explain how to obtain expressions for the higher order 
auto correlation functions of a density field from the mo- 
ments of the distribution of counts-in-cells. We also intro- 
duce the concept of higher order cross correlation functions, 
which combine fluctuations in two density fields. The con- 
cept of hierarchical amplitudes, scaling relations between 
higher order correlation functions and the two-point cor- 
relation function, is introduced in §2.2. The key theoretical 
results relating the higher order cross correlation functions 
of haloes to the two-point function and hierarchical ampli- 
tudes of the dark matter are given in §2.3. The simulations 
we use to measure the clustering of dark matter haloes are 
described in §2.4. 



2.1 The counts in cells approach to measuring 
clustering 

Here we give a brief overview of the approach of using the 
distribution of counts in cells to estimate the higher or- 
der auto correlation functions of a set of objects. An ex- 
cel lent and comprehensive review of this material is given 
bv lBernardeau et alj i|2002l ). We first discuss the higher or- 
der correlation functions for the case of a continuous, un- 
smoothed density field, then introduce the concept of cross- 
correlations (§2.1.1), before explaining how these results are 
changed in the case of a smoothed distribution of discrete 
points (§2.1.2). 



2.1.1 Higher order correlations: unsmoothed and 
continuous density field 

In general, the complete hierarchy of ./V-point correlation 
functions is required to fully characterize the spatial distri- 
bution of fluctuations in a density field. An exception to 
this occurs for the special case of a Gaussian density field, 
which can be described completely by its two-point correla- 
tion function. 

The iV-point correlation functions are usually written 
in terms of the dimensionless density fluctuation or density 
contrast at a point: 



5{x) = p(x)/(p) - 1, 



(1) 



where (p) is the mean density; the average is taken over 
different spatial locations. By definition, {8(x)) = when the 
average is taken over a fair sample of the density field. The 
7V th order moment of the density field, sometimes referred 
to as a central moment, because 8 is a fractional fluctuation 
around the mean density, is given by: 



pN 



(S(xi 



,8(Xn)), 



(2) 



where, in general, the density fluctuations are correlated at 
different spatial locations. 

The N th order central moments defined in Eq. [2] can 
be decomposed into terms which include products of lower 
order moments. This is because there are different permuta- 
tions of how the JV-points can be "connected" or joined to- 
gether. Th is idea is illustrated nicely by tree diagrams in the 
review by iBernardeau et all |2002h . The terms into which 
the central moments are broken down are called connected 
moments and these cannot be reduced further. In the tree 
diagram language, an TV-point connected moment has no 
disjoint points; all JV-points are linked to one another when 
the spatial averaging is performed. The distinction between 
connected and unconnected moments may become clearer if 
we write down the decomposition of the unconnected central 
moments up to fifth order: 



<<5 2 > = {8 2 ) C + (S) 2 C 

(5 3 ) = {8 3 } c + 3(5 2 } c {8) c + (8} s c 

(5 4 ) = {S 4 ) C + 4(S 3 ) C {S) C + 3(S 2 ) 2 C 

+6{S 2 ) c {5) 2 c + {S) i c 

(5 5 ) = {S 5 ) C + 5(S 4 ) C {S) C + W(S 3 ) C (S 2 ) C 

+10{5 3 ) C (5) 2 C + 15{5 2 ) 2 (8) C 
+W(5 2 ) C (8) 3 C + (S) 5 C , 



(3) 
(4) 
(5) 

(6) 
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where the subscript c outside the angular brackets denotes a 
connected moment. Remembering that (5) = 0, these equa- 
tions simplify to: 



Hence, for the second and third order moments, there is 
no difference in practice between the connected and uncon- 
nected moments. 

The JV-point auto-correlation functions, fjv, are written 
in terms of the connected moments: 



£n(xi, ■ ■ -,x N ) = (S(xi), . . .,S(x N )) c 



(11) 



By analogy with the TV-point auto correlation functions 
of fluctuations in a single density field, we can define the i+j- 
point cross correlation function of two, co-spatial density 
fields, with respective density contrasts given by 5i and 82: 



2/1, ■ --Vj) = 

(5i(xi), . . . , 5i(xi) 8 2 (yi) 



(12) 



><*2(j/j))c 



In the application in this paper, the first index will refer 
to the distribution of dark matter haloes and the second 
index to the dark matter. When the density contrasts are 
evaluated at the same spatial location, i.e. xi = . . . = Xi = 
yi — . . . — yj = 0, the connected moments £jj are called 
cumulants of the joint probability distribution function of 
Si and 82 (and are sometimes denoted as kij). 

To generate expressions for the higher order correlation 
functions of the cross-correlated density fluctuations, 
we will use the method of generating functions (see §3.3.3 
of Bernardeau et al. 2002). A moment generating function 
is defined for the central moments {fH,j) as a power series in 
81 and #2, which can be written as \ = (exp (5iii + (^2)}, 
where ti and t2 are random variables. This moment generat- 
ing function can be related to the cumulant generating func- 
tion (if)) for the connected cumulants by (see Bernardeau 
et al. 2002 for a proof): 



■t/j(t 1 ,t 2 ) = lnx(ti,t 2 ). 



(13) 



Then, by taking partial derivatives of ip and \ evaluated at 
ti = ti = 0, one can "generate" the cumulants and moments: 



fijj (0) — 



d 



i+i 



dt\ di? 2 

i+3 



d 



^|tl=*2=0 
X|tl=*2=0 



(sisi) 



(14) 
(15) 



Following this method we can obtain expressions for the 
cross-correlation cumulants up to order i + j = 5, grouping 
terms of the same order: 



k 2 ,o 

k3,0 

fe,i 



Mi,i 

M3,0 
M2,l 



(16) 
(17) 

(18) 
(19) 



^4,0 


— M4,0 


n 2 
— 3/X2.0 


^3,1 


= M3,l 


— 3 /i2,oMi,i 


k2,2 


= M2,2 


2 
— ^2,0^0,2 — ^ 



(8 2 ) 


= (S 2 )c 


(7) 






(8 s ) 


= (5 3 ) c 


(8) 


ks,o 


= M5,0 




= {8 4 } C + 3(8 2 ) 2 C 


(9) 


k4,i 


= M4,l 


(8 5 ) 


= {8 3 ) C + W(S 3 ) C (S 2 ) C . 


(10) 


&3,2 


= M3,2 



(20) 
(21) 
(22) 



10 /X 3 , 0^2,0 (23) 

4 /U3,oA*i,i — 6 fi 2 ,ofJ-2,i (24) 

M3,0M0,2 - 6/X 2 , — 3^2,0Ml,2- (25) 

Note that these results are symmetric with respect to ex- 
changing the indexes and that we have used the fact that 
fJ.1,0 = /tio.i = 0, since, by construction (81) = (82) = 0. 



2.1.2 Higher order correlations: smoothed and discrete 
density fields 

Sadly, density fluctuations at a point are of little practi- 
cal use as they cannot be measured reliably, as typically 
we have a finite number of tracers of the density field, i.e. 
galaxies in a survey or dark matter particles in an N-body 
simulation, and so have a limited resolution view of the den- 
sity field. Furthermore, estimating the iV-point correlations 
for a modern survey or simulation is time consuming and 
short cuts are often taken, such as restricting the number 
of configurations of points sampled. To overcome both of 
these problems, moments of the smoothed density field can 
be computed instead of the point moments. 

The smoothed density contrast, Sr, is a convolution of 
the density contrast at a point with the smoothing window, 
Wr, which has volume V: 



S(x) R = -j^ / da; 3 8{x)Wr{x — 2 



(26) 



Typically, the smoothing window is a spherical top-hat in 
which case Wr = 1 for all points within distance R from the 
centre of the window and Wr — otherwise. After smooth- 
ing, the cumulants correspond to the i + j-point volume- 
averaged cross correlation functions: 



d 3 xi . 



.d 3 Xi 



dV 



Vj 



(27) 



W R {xx) . . . Wr{ Xi ) W R (vi) . . . W R ( yj )ii. 



Eqs. I16H25I are still valid, with the cumulants replaced by 
volume-averaged cumulants. 

Another issue introduced by the discreteness of the den- 
sity field is the contribution of Poisson noise to the measure- 
ments of the cumulants. To take this into account, we can 
modi fy the moment generating function as follows (|Peebled 
I1980T I: 



X (ti,t 2 ) = (exp(./i (fi) + / a (t 2 ))>, (28) 
h = (exp(ti) - ii - 1) ni + (exp(ti) - 1) <5i, (29) 
h = (exp(t 2 ) - *2 - 1) fta + (exp(t 2 ) - 1) £2. (30) 

Here, fi\ and fi2 are the mean number of objects in density 
field 1 and density field 2 respectively within spheres of ra- 
dius R. Using this modified generating function, and defining 
thj — (( n i — ^i)* (™2 — n 2 ) 3 ), we obtain the following rela- 
tions between the volume-averaged, connected i + j-point 
cross correlation functions, and the central moments, 



4 Angulo et al. 



n?6,o = /4,o-"i (31) 

nin 2 fi,i = n'i t i (32) 

n\ 6,2 = Mo, 2 ~ "-2 (33) 

n\ 6,o = M3,o + 2ni - 3m 2 ,o (34) 

n?n 2 6,i = M2,i - Mi,i (35) 

mn|fi,2 = m'i, 2 - (36) 

"•2 6,3 = Mo, 3 + 2n 2 - 3mo,2 (37) 

«i 6,0 = M4,o - 6ni + 11^2,0 - 6m 3 ,o - 3/x^o (38) 

n\n 2 l-i,i = M3,i + 2 m'i,i ^ 3 M 2 ,i - 3/ii,iM2,o (39) 

-2— 2 ^ / / / / ll l ac\\ 

nin 2 6,2 = M2,2 - Mi,2 - M2,i + Mi,i - M2,oMo,2 (40) 

o '2 

-2mi,i 

Win2 6,3 = Mi,3 + 2m'i,i - 3m'i,2 - 3/ii,iMo,2 (41) 

ni 6,4 = Mo,4 - 6n 2 + ll/4), 2 - 6^0,3 - 3mo% (42) 



Note that these expressions revert to those in the literature 
for autocorrelation mom ents in the case of e ither i or j equal 
to zero (see for instance iBaueh et al.|[l995h . Also note that 
in the limit fii — > 00, fi2 — » 00, they correspond to the 
expressions given by Eqs. 1 161251 

2.2 Hierarchical amplitudes 

At this point it is useful to define quantities called hier- 
archical amplitudes which are the ratio between the N- 
point, volume-averaged connected moments and the two- 
point volume-averaged connected moment raised to the 
TV — 1 power: 

Sn = J^. (43) 

?2 

This form is motivated by the expected properties of a 
Gaussian field which evolves due to gravitational instabil- 
ity (Bernardeau et al. 2002). In the case of small amplitude 
fluctuations, i.e. on smoothing scales for which fe(R) <C 1, 
the Sn depend only on the local slope of the linear pertur- 
bation theory power spe ctrum of density fluctuations and 
are in depe ndent of time Jjuszkiewicz, Bouchet. fc Colombil 
1 19931 ; see iBernardeaul 1 19941 for expressions for the Sn). 
Similar scalings, but with different values for the Sn, ap- 
ply in the case of distributions of particles which have 
not arisen through gravitational instability, e.g. parti- 
cles displaced according to Zeldovich approximation (see 
lJuszkiewicz. Bouchet. fc Colombilll993l ). 

In the case of a Gaussian density field, all of the Sn 
are equal to zero. Initially, as perturbations grow through 
gravitational instability, the two-point connected moment 
increases. The distribution of fluctuations soon starts to de- 
viate from a Gaussian, particularly as voids grow in size 
and cells become empty (S — > —1). Voids evolve more slowly 
than overdense regions. There is in principle no limit on 
how overdense a cell can become. As a result, the distribu- 
tion of overdensities becomes asymmetrical or skewed, with 
the peak of the distribution moving to negative density con- 
trasts and a long tail developing to high density contrasts. 
To first order, this deviation from symmetry is quantified by 



the value of S3, which is often referred to as the skewness 
of the density field. Higher order moments and hierarchical 
amplitudes probe progressively further out into the tails of 
the distribution of density contrasts. 

2.3 Higher order correlations: biased tracers 

We are now in a position to consider the cross-correlation 
functions for the case of relevance in this paper, when the 
set of objects making up one of the density fields is local 
function of the second density field; the first density field is 
a biased tracer of the second. In our application, one density 
field is defined by the spatial distribution of dark matter 
haloes and the other by the dark matter. In the case of a 
local bias and small perturbations, the density contrast in 
the biased tracers (Si) can be written as an expansion in 
terms of the underlying dark matter density contrast (82), 
as proposed by Fry & Gaztanaga (1993): 
00 

fc=0 

where the bk are known as bias coefficients; bi is the linear 
bias commonly discussed in relation to two point correla- 
tions. Note that, by construction, we require that (S) = 0, 
which implies bo = — '^2^- 2 (pk) /k\. The bk, as we shall see 
later, depend on mass but this is suppressed in our notation. 

Using this bias pr e script ion, and following the treat- 
ment iFrv fc Gaztanaeal j 1993ft used for autocorrelations, we 
can write the volume-averaged cross-correlation functions of 
dark matter haloes in terms of the two-point volume aver- 
aged correlation function (6,2) and hierarchical amplitudes 
of the dark matter, Sn: 



6,1 = &i£o,2 + O (g, a ) (45) 

6.0 = &?fo,2 + O (g, a ) (46) 

£i,a = 6ig,a (c 2 + S 3 ) + O (g, a ) (47) 

6.1 = &?S,2 (2 c 2 + S3) + 0(^,2) (48) 

6.0 = &iCL(3c2 + S3)+0(fo, 2 ) (49) 

6.3 = htn, 2 (3S 3 c 2 + S 4 + c 3 ) + O (Co, 2 ) (50) 

6.2 = &?fo,2 (Si + 6S 3 c 2 + 2c| + 2c 3 ) + O (6 4 , 2 ) (51) 

6.1 = &?S>,2 (6(1 + 9S 3 c 2 + S 4 + 3c 3 ) + O (6 4 , 2 ) (52) 
6,0 = Bj*o,2 (I2ca + 12S 3 c 2 + S 4 + 4c 3 ) + O (£0,2) (53) 

6.4 = bi6 4 ,2(4c2S4 + 6c 3 S 3 + (54) 



c 4 + S 5 + 3c 2 Sf ) + O (6 5 ,2) 

6,3 = &?f£ 2 (12S 3 C 3 + 6Sf C 2 + 12S 3 c| + 6c 2 C 3 + (55) 

2c 4 + S 5 + 8c 2 S 4 ) + O (g) 

6,2 = &?6 4 ,2(12C2S4 + 18C3S3 + 18C2C3+36c|S 3 (56) 

+9c 2 S 3 2 + S 5 + 6ci + 3c 4 ) + O (6 5 ) 
6,1 = &ifo,s(4c4 + 24ci + S5 + 72clS3 + 16c2S 4 + (57) 
36c 2 c 3 + 24c 3 S 3 + I2C2S3) + O (6 s ) 
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6,0 = &i^o,2(20c25 , 4 + 15c 2 5'3 2 + 60ci + 30c3S3+ (58) 

5c 4 + 120c|5* 3 + Ss + 60c 2 c 3 ) + O 

where Ck = bk/bi. Note it has been shown that these trans- 
formations preserve the hierarchical nature of the clustering 
|Frv fc Gaztanagalll993| y 



2.4 Numerical Simulations 

To make accurate measurements of the higher order clus- 
tering of dark matter and dark matt er haloes, we use th e 
N-body simulations carried out by lAngulo et alj |2008). 
Two simulation specifications were used: i) The BASICC, 
a high-resolution run which used 1448 3 particles of mass 
5.49 x 10 11 h~ x Mq to follow the growth of structure in the 
dark matter in a periodic box of side 1340/i _1 Mpc. ii) The 
L-BASICC ensemble, a suite of 50 lower resolution runs, 
which used 448 3 particles of mass 1.85 x 10 12 h' 1 M Q in the 
same box size as the BASICC. Each L-BASICC run was 
evolved from a different realization of the initial Gaussian 
density field. The simulation volume was chosen to allow the 
growth of fluctuations to be modelled accurately on a wide 
range of scales, including that of the baryonic acoustic oscil- 
lations (the BASICC acronym stands for Baryonic Acoustic 
oscillation Simulations at the Institute for Computational 
Cosmology). The extremely large volume of each box also 
makes it possible to extract accurate measurements of the 
clustering of massive haloes. The superior mass resolution 
of the BASICC run means that it can resolve the haloes 
which are predicted to host the galaxies expected to be seen 
in forthcoming galaxy surveys. The L-BASICC runs resolve 
haloes equivalent to group-sized systems. The independence 
of the L-BASICC ensemble runs makes them ideally suited 
to the assessment of the impact of cosmic variance on our 
clustering measurements. 

In both cases, the same values of the basic cosmological 
parameters were adopted, which are broadly consistent with 
recent data from the cosmic microwav e background and th e 
power spectrum of galaxy clustering (jSanchez et al.l i2006): 
the matter density parameter, Qm = 0.25, the vacuum en- 
ergy density parameter, £7a = 0.75, the normalization of 
density fluctuations, expressed in terms of the linear the- 
ory amplitude of density fluctuations in spheres of radius 
8/i _1 Mpc at the present day, as = 0.9, the primordial spec- 
tral index n s = 1, the dark energy equation of state, w = — 1, 
and the Hubble constant, h — //(>/( 100kms _1 Mpc _1 ) = 
0.73. The simulations were started from realizations of a 
Gaussia n density field se t up using the Ze'ldovich approxi- 
mation ()Zel Particles were pertur bed from a 
glass-like distribution (|Baugh et aTll 19951 : [Whit3l 1994^ . The 
starting redshift for both sets of simulations was z — 63. 
The linear perturbation theory power spectrum used to set 
up the initi al density field was generated using the Boltzman 
code CAMB ()Lewis et al.l l2000). The initial density field was 
evolved to th e present day using a memory efficient version 
of GADGET-2 (|Springelll2005h . 

Outputs of the particle positions and velocities were 
stored from the simulations at selected redshifts. Dark mat- 
ter haloes were identi fied using a Friend s-of-Friends (FOF) 
percolation algorithm ijDavis et al.ll 19851 ) and substructures 
within these were fou nd using a modified version of SUBFIND 
l|Springel et al.|[200ll ). Our default choice is to use the num- 
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Figure 1. The ratio V ma x/V200 as a function of halo mass for 
gravitationally bound haloes in the BASICC simulation, which have 
a minimum of 26 particles. V m ax is the maximum effective circular 
velocity of the largest substructure within the halo and V200 is 
the effective rotation speed at the radius within which the mean 
density is 200 times the critical density, computed using all of 
the particles within this radius. Each panel shows the relation 
at a different redshift as indicated by the legend. The red lines 
show the 20-80 percentile range of the distribution of Vmax/Vioo 
values, and the blue lines show the mean. 



ber of particles in a structure as returned by the FOF group 
finder to set the mass of the halo; at the end of Section 3.4 
we discuss a variation on this to assess the sensitivity of 
our results to the group finder. The position of the halo is 
the position of the most bound particle in the largest sub- 
structure, as determined by SUBFIND. In this paper, only 
gravitationally bound groups with more than 26 particles 
are considered. The SUBFIND algorithm also computes sev- 
eral halo properties such as the circular velocity profile 
V c (r) = (GM(r)/r)^ 2 , V max , the maximum value of V c for 
the largest substructure, and V^oo = Vc(r20o), where r2oo is 
the radius of a sphere enclosing a volume of mean density 200 
times the critical density. These properties are calculated us- 
ing only the particles which are bound to the main subhalo 
of the FOF halo; i.e. ignoring all of the other substructure 
haloes within the FOF halo. In the best resolved haloes, sub- 
structures other than the largest substructure accoun t for at 
most 15% of the total halo mass l|Ghigna et al.|[l99ct ). Later 
on in the paper we will present results for the clustering of 
haloes as a function of mass and a second parameter. We 
have a limited number of output times available to us, so it 
is not feasible to use the formation time of the halo as the 
second parameter. Instead, we will use the ratio Vmax/V^oo- 
Fig. 1 shows Vmax/V^oo as a function of halo mass at differ- 
ent epochs in the BASICC simulation. There is a trend of de- 
clining Vmax/V^oo with increasing halo mass. In cases where 
the density profile of the d ark matter halo matc hes the uni- 
versal profile advocated bv lNavarro et all l|l997l ). Knax/V^oo 
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depends on the concentration parameter which characterizes 
the profile. Haloes in the extreme parts of the distribution 
of Vmax/V20 also have extreme values of th e concentration 
parameter (|Navarro. Frenk. fe Whitel Il997l ) . More massive 
haloes tend to have lower values of the concentration pa- 
rameter and lower values of the velocity ratio V ma x/V^oo- 
The ratio Vmax/Vwo is easier to extract from the simula- 
tion, as it does not require a parametric form to be fitted 
to the density profile. There is a correlation between forma- 
tion time and concentration para meter, and hence the ra- 
tio Vmax/V200i albeit w ith scatter (|Navarro. Frenk. fc White! 
ll997l ; IZhao et al.ll2003l ). 



3 RESULTS 

Our ultimate goal is to measure the higher order bias of 
dark matter haloes. As described in Section 2, we follow 
a novel approach to do this, employing cross moments be- 
tween haloes and the dark matter. The first step in this pro- 
cess is to compute the densities of haloes and dark matter 
on grids of cubical cells of different sizefl A natural by- 
product of this procedure is the higher order clustering of 
the dark matter and haloes in terms of the auto-correlation 
functions. We first present the hierarchical amplitudes esti- 
mated for the dark matter (§3.1) and haloes (§3.2) using the 
autocorrelation function higher order moments. In §3.3 we 
show the measurements of the cross moments and in §3.4 
we present the interpretation of these results in terms of the 
bias parameters. 

3.1 Hierarchical amplitudes for the dark matter 

Fig. [5] shows the hierarchical amplitudes Sn measured for 
the dark matter at different redshifts. The upper panels 
show the results in real-space and the lower panels in- 
clude the effects of redshift space distortions using the dis- 
tant observer approximation. The points indicate the me- 
dian value of the hierarchical amplitudes measured in the 
L-BASICC ensemble and the error bars indicate the vari- 
ance in these measurements. The lines show the hierarchical 
ampl i tudes predicted by perturbation theory (iBernardeaul 
1 19941 ; lJuszkiewicz. Bouchet. fc Colombill993 ). At the high- 
est redshift plotted, 2 = 4, the agreement between the mea- 
surements made from the simulations and the predictions 
of perturbation theory is impressive, covering scales from 
5h _1 Mpc to 100/i _1 Mpc for S 3 and S 4 . As redshift de- 
creases, the simulation results for S5 and 56 are slightly 
higher than the perturbation theory predictions. The mea- 
surements of S3 from the simulations continue to agree with 
the perturbation theory predictions, but over a narrower 
range of scales. For smoothing scales on which the variance is 
less than unity, the hierarchical amplitudes are expected to 
be independent of epoch, depending only on the shape of the 
linear perturbation theory power spectrum of d e nsity fluctu- 
ations (jjuszkiewicz. Bouchet. fc Colomb i 1993; B ernardeaul 

1 Tests show that density fluctuations in cubical cells can be read- 
ily translated into counts in spherical cells by simply setting the 
volume of the spherical cell equal to that of the cube. We use 
cubical cells for speed. The counts are regridded to improve the 
measurement of the rare event tails of the count distribution. 



1 19941 ; iGaztanaga fc Baugh|[l995l 'l. Fig. [2] confirms that this 
is the case. As the density field evolves, the measured hi- 
erarchical amplitudes change remarkably little, particularly 
when one bears in mind that the higher order correlation 
functions change substantially between z = 4 and z = 0. 
For example, for a cell of radius 50/i -1 Mpc, the two-point 
volume averaged correlation function increases by a factor 
of 14 over this redshift interval, and the three-point function 
by a factor of 197. 

Nevertheless, the simulation results do tend to exceed 
the perturbation theory predictions on all scales at all orders 
as the density fluctuations grow. 

The hierarchical amplitudes measured on small scales 
differ significantly from the predictions of perturbation the- 
ory. At z = 4, the simulation results are below the analyt- 
ical predictions for cell radii smaller than R ~ 5/i _1 Mpc. 
This behaviour is sensitive to the arrangement of particles 
which is perturbed to set up the initial density field. At later 
times, the memory of the initial conditions is erased on small 
scales and the measured amplitudes greatly exceed the ex- 
pectations of perturbation theory. On these scales, the dom- 
inant contribution to the cross correlation moments is from 
particles within common dark matter haloes. Note that in 
Fig. [2] we do not correct the measured higher order correla- 
tion functions for Poisson noise, since the initial density field 
was created by perturbing particles distributed in a glass- 
like configuration which is sub-Poissonian. Hence, the dark 
matter density field is not a random samp ling of a contin- 
uous density field fsee lAngulo et all l|2008l ) for an extended 
discussion of this point). The turnover in the hierarchical 
amplitudes seen at small cell radii (e.g. for R < 2/i _1 Mpc is 
due to the finite resolution of the L-BASICC simulations; the 
hierarchical amplitudes continue to increase in amplitude on 
smaller smoothing scales in the BASICC run. 

The lower panels of Fig. [5] show the impact of gravi- 
tationally induced peculiar motions on the hierarchical am- 
plitudes. We model redshift space distortions using the dis- 
tant observer approximation, in which peculiar motions per- 
turb the particle position parallel to one of the co-ordinate 
axes. Virialized structures appear elongated when viewed in 
redshift space. On large scales, coherent bulk flows tend to 
increase the amplitude of correlation functions. There is a 
modest reduction in the amplitude of the hierarchical am- 
plitudes on large scales. On small scales, there is a dramatic 
reduction in the magnitude of the Sn. The overall impact 
of the redshift space distortions is to greatly reduce the de- 
pen dence of the hierar chical amplitudes on smoothing scale 
(see iHovle et al.1 12000|) . 

The estimated error on the measured hierarchical mo- 
ments is shown in Fig. [3] in which we plot the fractional 
error on Sn obtained from the scatter in the measurements 
from the L-BASICC ensemble. The plot suggests that the 
skewness of the dark matter can be well measured on all 
smoothing scales considered from a volume of the size of the 
L-BASICC simulation cube. The range of scales over which 
robust measurements can be made of the hierarchical ampli- 
tudes becomes progressively narrower with increasing order. 
For example, at z = 0, reliable measurements of Se are lim- 
ited to smoothing radii smaller than R ~ 30/i _1 Mpc. 
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Figure 2. The hierarchical amplitudes (Sn) measured for the dark matter as a function smoothing scale, which is plotted in terms of 
the radius of the sphere with the same volume as the cubical cell used. The upper panels show the results in real-space and the lower 
panels show redshift-space. Each panel corresponds to a different redshift as indicated by the legend. The points show the amplitude for 
the Sn obtained from the L-BASICC ensemble, after taking the ratio of the median correlation functions, as defined by Eq. I43I The error 
bars show the scatter in the measurements over the ensemble, obtained by computing Sn for each simulation from the ensemble. Error 
bars are plotted at smoothing scales for which the fractional error is less than unity; triangles show scales on which the error exceeds 
unity. In both sets of panels, the dashed lines show the predictions of perturbation theory in real space (see text for details). Note that 
no correction for shot noise has been applied to the measured amplitudes. The arrows indicate the cell radius for which the variance in 
the counts in cells for the dark matter is equal to unity, which is roughly the scale down to which perturbation theory should be valid; 
at z = 4. this scale is below R = 1/i _1 Mdc. 
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Figure 3. The fractional scatter, <r(Ss)/Stq, in the measured hierarchical amplitudes, as estimated from the 50 simulations in the 
L-BASICC ensemble. Different lines show the scatter for different orders as indicated by the legend. 



3.2 The hierarchical amplitudes of dark matter 
haloes 

The hierarchical amplitudes of dark matter haloes are more 
complicated than those of the dark matter. In addition to 
a term arising from the evolution of the density field un- 
der gravitational instability, there is a contribution which 
depends upon height of the peak in the initial density fiel d 
which collapses to form the halo l|Mo. Jing. fc Whitell 19971 ). 
For example, if we consider the second and third order auto- 
correlation functions of haloes given by Eqs. 44 and 47, then 
the skewness for dark matter haloes, S^ , is given by: 



sf 



£3,0 

(&.o) a 
362 Ss 
b\ + 61 ' 



(59) 
(60) 



The gravitational contribution to the skewness, S3, is di- 
luted by the linear bias factor, 61 . In the case of rare peaks, 
or, equivalently, haloes with masses far in excess of the char- 
acteristic mass, Al t , at a given redshift, approaches an 
asymptotic value. In this limit, ~ and so S3 1 ~ 3; simi- 
lar arguments for the fourth and fifth order hierarchical am- 
plitudes yield asymptotic values of = 16 and Sf — 125 
(|Mo, Jing. fc White! fl997h . Massive haloes at high redshift 
can therefore have non-zero hierarchical amplitudes even if 
the dark matter distribution still has a Gaussian distribution 
and hence 5° M = 0. 

We plot the hierarchical amplitudes of dark mat- 
ter haloes in Fig. [4] as a function of the scaled peak 



height, S c /a(M, z). The simulation results are averaged over 
smoothing radii of 20fr~ 1 Mpc < R < 50/i _1 Mpc. The 
dashed line sho ws the prediction obtaine d assuming the 
mass function of iPress fc Schechterl (Il974l) and the spher- 
ical collapse model (see lMo. Jing, fc Whitell 19971 ) . The solid 



line shows an improved calculation which uses t he ellipsoidal 
collap se model and the mass function derived bv lSheth et al.l 
|200ll ). There is some dispersion between the simulation re- 
sults at different redshifts. The measurements are in rea- 
sonable agreement with the theoretical predictions for large 
values of S c /a(M, z). For more modest peaks, the hierarchi- 
cal amplitudes of haloes averaged on large smoothing scales 
show a dip and are significantly smaller than the amplitude 
recovered for the dark matter. The strength of this dip is 
more pronounced in the measurements from the simulations 
than it is in the theoretical predictions. This discrepancy 
suggests that the theoretical models do not reproduce the 
trend of bias with halo mass for such objects, as we shall see 
in §3.4. 



3.3 Cross-correlation estimates of higher order 
clustering 

We now switch to estimating cross-correlation functions in- 
stead of auto-correlation functions. To recap §2, to reduce 
the impact of discreteness noise on our measurement of halo 
clustering, we cross-correlate fluctuations in the spatial dis- 
tribution of haloes with the fluctuation in the dark mat- 
ter density within the same cell. As the order of the corre- 
lation function increases, the number of possible permuta- 



The assembly 



bias of dark matter haloes to higher orders 9 



300 
250 
200 
150 
100 
50 


-50 



: □ z = 0.0 

L X z = 0.5 
: Oz = 1.0 
: A z = 2.0 

'- X z = 3.0 

>• > 












- B — □ • □ 





20 






15 


: > 




10 








L 







" - □ _Q- 




-5 






-10 








<5 c /<7(M,z) 

Figure 4. The hierarchical amplitudes of dark matter haloes, 
plotted as a function of the peak height corresponding to the 
halo mass. In this plot, the hierarchical amplitudes are averaged 
over cell sizes of 20?t -1 Mpc < R < 50/i -1 Mpc. The dashed curve 
shows a theoretical prediction based on the spherical collapse 
model (Mo, Jing & White 1997) and the solid line shows a re- 
vised prediction based upon an ellipsoidal collapse, by Sheth, Mo 
& Tormen (2001). The corresponding hierarchical amplitudes for 
the dark matter, averaged over the same range of cell radii, are 
indicated in each panel by the arrow. 



tions of halo fluctuations and dark matter fluctuations in- 
creases. For a given order of correlation function, the rela- 
tion between these permutations can be understood using 
the expressions for the cross-moments given in §2.3. The 
relationship at second order is particularly straightforward. 
The halo autocorrelation function, £2,0 (recall the first index 
gives the order of the halo density contrast and the second 
index gives the order of the dark matter density contrast) is 
related to the auto-correlation of the dark matter, £0,2, by 
£2,0 = &i£o,2- The second-order cross correlation function, 
£1,1, is related to the autocorrelation function of dark matter 
by £1,1 = &i£o,2- The primary difference between £2,0 and £1,1 
is therefore a factor of 61 . This basic trend is approximately 
replicated for any order of correlation function: as fluctua- 
tions in the halo density are substituted by fluctuations in 
the dark matter, the amplitude of the cross-correlation is 
reduced by a factor which depends on 61 . Above second or- 
der, this factor is modulated by higher order bias terms and 
the hierarchical amplitudes of the dark matter (see §2.2). 
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Figure 5. Volume-averaged i + j-point cross-correlation 
functions, measured for haloes of mass 1.1 X 10 13 < 

(Al/h- 1 M Q ) < 2.8 X 10 13 (619386 objects) and the dark matter 
at 2 = in the BASICC simulation. The auto-correlation func- 
tion of haloes is denoted £i+j,o and the auto-correlation of dark 
matter by £o,i+j- Each panel shows a different order of cross- 
correlation. The key shows the different permutations of cross- 
correlation function in each case. The moments have been cor- 
rected for Poisson noise due to the finite number of haloes. 
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Figure 6. The fractional error on the four-point cross correla- 
tion functions, estimated from the scatter over the L-BASICC runs. 
Each panel shows the results for a different redshift, as shown by 
the key. The legend shows the different permutations of cross- 
correlation moment. To improve the statistics, all the haloes in 
the L-BASICC runs have been used in this case. 



The precise relation between the different permutations of 
cross-correlation functions depends upon the values of the 
bias parameters and therefore on the halo mass under con- 
sideration. 

We show illustrative examples of volume-averaged 
cross-correlation functions, £i t j, estimated from the BASICC 
simulation in Fig. [S] Each panel shows a different order 
of clustering, starting with the second moment in the top 
left panel and ending with the fifth order correlation func- 
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tion in the bottom right panel. In this plot, the haloes used 
have masses in the range 1.1 x 10 13 < (Mhaio/ft. -1 Mq) < 
2.8 x 10 13 and the clustering is measured at z = 0. The 
top-left panel of Fig. [5] shows that there is little difference 
in the amplitude of the second-order correlation function on 
large smoothing scales between the different permutations 
of This implies that for these haloes, the linear bias 
term bi w 1. The correlation functions are, however, differ- 
ent on small scales. The autocorrelation function of the dark 
matter (£0,2) is steeper than the autocorrelation of haloes 
(£0,2). The cross-correlation functions are different on large 
scales for third, fourth and fifth orders. The difference in 
amplitude is fairly independent of scale for cells with radii 
R > 10/i _1 Mpc. Since the linear bias of this sample of haloes 
is close to unity, this difference is driven by the higher or- 
der bias terms and the hierarchical amplitudes of the dark 
matter. We plan to model the full behaviour of the cross- 
correlation functions, including the small scale form, using 
the halo model in a future paper. 

One might be concerned that replacing fluctuations in 
halo density by fluctuations in dark matter in the higher 
order correlation functions leads to a reduction in the clus- 
tering amplitude (as is indeed apparent in Fig. [5]). How- 
ever, this is more than offset by a reduction in the noise 
or scatter of the measurement. The fractional error on the 
measurements of the cross-correlation functions is plotted 
in Fig. [6] The scatter is estimated using the L-BASICC en- 
semble. Each panel shows the scatter at a different redshift. 
The cross correlation (i.e. one part halo fluctuation, 

i + j — 1 parts dark matter fluctuation) gives the optimal 
error estimate, with a performance comparable to the auto- 
correlation of the dark matter. At z = 1, it is not possible to 
measure the four-point autocorrelation function of this sam- 
ple of haloes, even with a box of the size of the L-BASICC 
runs. Nevertheless, it is possible to measure the bias factors 
relating the four-point functions of haloes and mass using 
the cross-correlation. Our use of a cross-correlation estima- 
tor therefore allows us to extend the measurements of the 
higher order clustering of haloes to orders and redshifts that 
would not be possible using auto-correlations. 



3.4 The bias parameters of dark matter haloes 

We now use the cross correlation functions to estimate the 
linear and higher order bias parameters of dark matter 
haloes. As we demonstrated in the previous section, the best 
possible measurement of the i + j order correlation func- 
tion is obtained when the cross correlation function is made 
up of one part fluctuation in halo density and i + j — 1 parts 
dark matter fluctuation: i.e. in our notation This 
approach, combined with the huge volume of our simula- 
tion, makes it possible, for the first time, to measure the 
third and fourth order bias parameters, and to do so using 
narrow mass bins. 

In this section, we use the higher resolution BASICC run, 
which can resolve the largest dynamic range in halo mass. 
We use the higher order correlation function measurements 
over the range of smoothing radii 15 < (R/h~ 1 Mpc) < 50 
to estimate the halo bias parameters. The large volume of 
the BASICC simulation means that we can make robust mea- 
surements of the higher order correlation functions out to 
larger smoothing radii than is possible with the smaller Mil- 



lennium simulation. The smallest scale we use is set by the 
requirement that the expansion relating the overdensity in 
haloes to the overdensity in dark matter (Eq. 44) is a good 
approximation, i.e. when £ <C 1. The scales we use to extract 
the halo bias parameters are considerably larger than those 
Gao, Springel & White (2005) and Gao & White (2007) 
were able to use in the Millennium. We use the simulation 
outputs at redshifts of z = 0, 0.5, 1, 2 and 3 to measure the 
clustering of haloes. 

The results for the first, second, third and fourth or- 
der bias parameters of dark matter haloes are presented in 
Fig. [7] Each panel corresponds to a different order. The up- 
per half of each panel shows the respective order of bias 
parameter as a function of halo mass, expressed in terms of 
the peak height corresponding to the halo mass, S c /a(M, z). 
The lower half of each panel shows the deviation from the 
bias parameter extracted for a given mass for samples of the 
20% of haloes in the mass bin with the highest and low- 
est values of Vmax/V^oo, which we are using as a proxy for 
halo concentration. Different symbols in the upper panels 
show the measurements at different output redshifts in the 
BASICC run, as indicated by the key; the same colours are 
used to draw the lines showing results for samples defined 
by different Kn ax /V2oo values at the same output redshifts 
in the lower panels. 

In Fig. [7] there is remarkably little scatter between the 
results obtained from the different output redshifts for the 
case of the overall bias as a function of mass. This is en- 
couraging, as it shows that our results are not affected by 
resolution (haloes with similar values of 5 c /a(Al, z) at dif- 
ferent output times are made up of very different numbers of 
particles). Gao, Springel & White (2005) were able to mea- 
sure the linear bias parameter up to haloes corresponding 
to peak heights of 3<r; we are able to extract measurements 
for haloes corresponding to 5cr peaks. 

In the upper sub-panels of Fig. [7] we show two the- 
oretical predictions for the bias parameters of dark mat- 
ter haloes. The dotted lines show the predictions from 
iMo. Jing. fc White! l|l997h . based on an extension of Press 
& Schechter's (1974) theory for abundance of dark mat- 
ter haloes and the spheric al collapse model. The sol id lines 
show the calculation from Scoccimarro ct al. (2001) which 
uses the mass function of ISheth fc Tormenl (|l999h . Our re- 
sults tend to agree best with the later, although the mea- 
surements favour a steeper dependence of bias on peak 
height at all orders. For less rare peaks, neither theoret- 
ical model gives a particularly good fit to the simula- 
tion results. A similar trend, albeit with more scatter be- 
tw een the results at differ e nt ou tput reds hifts was found 
by iGao. Springel. fc White! l[2005l) (see also IWechsler et~aH 
120061 : iJing. Suto. fc Moll2007n . 

Previous studies have reported a dependence of cluster- 
ing strength on a second halo property besides mass, such as 
halo formation time or concentration (|Wechsler et al. I l2006l ; 
Gao & White 2007). We do not have sufficient output times 
to make a robust estimate of formation time so we use a 
proxy for halo concentration instead, Knax/Vajo- We find 
that the clustering of high peak haloes is sensitive to whether 
the halo has a high or low value of Vmax/V^oo- The 20% of 
haloes with the lowest values of V^ax/V^oo within a given 
mass bin (i.e. those with the lowest concentrations) have 
the largest linear and second order bias terms. This result 
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Figure 7. The bias parameters as a function of halo mass parametrized by v = 5 c /a(M, z). Each plot shows a different order of bias 
parameter: a) linear bias 61, b) the ratio of the 2nd order bias, 62/61, c) the ratio of the 3rd order bias, 63/61 and d) 64/61. In the 
lower panel of each plot, the residual bias parameters for the 20% of haloes with the highest or lowest values of V max /V200i a proxy for 
concentration, are plotted. In the upper panels, symbols show the measurements for different output redshifts, as indicated by the key. 
The same line colours are used to show the results for different redshifts in the lower panels. In the upper panel of each plot, we plot two 
theoretical predictions for the bias parameters, given by Mo, Jing & White (1997) and Scoccimarro et al. (2001). 



agrees with previous estimates of the dependence of the lin - 
ear bias term on halo concentration ()Wechsler et al-l feoOG), 

The peak height dependence of the third and fourth 
order bias terms for haloes split by V ma x/V2oo is more com- 
plicated. Fig. [7] shows that the third order bias depends on 
our concentration proxy in a non-monotonic fashion. The 
trend for the fourth order bias is reversed compared with 
the results for the first and second order bias parameters: 



low concentration haloes have a negative value of the fourth 
order bias. We note that it would not be possible to mea- 
sure a fourth order bias at all using halo auto-correlation 
functions. 

One might be concerned that our results could be sen- 
sitive to the operation of the group finder. In particular, it 
is well known that the FoF algorithm can sometimes spuri- 
ously link together distinct haloes into a larger halo, through 
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bridges of particles (e.g. ICole fc Lacevlll996l ). We therefore 
carried out the exercise of relabelling the mass of each halo 
by the mass of the largest substructure as determined by 
SUBFIND. In the rare cases in which haloes are incorrectly 
linked into a larger structure, using instead the SUBFIND 
mass would result in a significant shift in the mass bin to 
which the halo is assigned. Moreover, one would expect that 
low concentration FoF haloes would be more prone to being 
broken up in this way. However, we found no change in our 
results upon following this procedure, demonstrating that 
the trends we find for the dependence of bias on mass and 
concentration are robust. 



4 SUMMARY AND DISCUSSION 

In this paper, we have combined ultra-large volume cosmo- 
logical simulations with a novel approach to estimating the 
higher order correlation functions of dilute samples of ob- 
jects. The large simulation volume allows us to extract bias 
parameters on large scales, which follow linear perturba- 
tion theory more closely, and provides us with large samples 
of high mass haloes from which robust clustering measure- 
ments can be made. The cross-moment counts-in-cells tech- 
nique we use to estimate the higher order clustering of dark 
matter haloes has superior noise performance to traditional 
autocorrelation functions, allowing us to probe clustering 
to higher orders. These improvements made it possible to 
extend previous work on the assembly bias of dark matter 
haloes in a number of ways. We have been able to extract 
measurements of halo clustering for objects corresponding to 
5<t peaks, almost twice as high as in earlier studies. We have 
also presented, for the first time, estimates of the higher or- 
der bias parameters of haloes, up to fourth order, and using 
narrow mass bins. 

Our results are in qualitative agreement with those in 
the literature where they overlap. We find that the linear 
bias factor, bi, is a strong function of mass, varying by 
an order of magnitude for peaks ranging in height from 
S c /a(M, z) — 1 to 5. We use the ratio of the maximum 
of the effective halo rotation speed to the speed at the virial 
radius, V ma x/V2oo as a proxy for halo concentration. High 
mass, high Vmax/Vajo haloes are less strongly clustered than 
the same mass haloes with low values of Knax/Vioo; haloes 
with S c /a(M, z) ~ 4 display second order clustering that 
differs by ~ 25% between the 20% with the lowest values of 
Vinax/V^oo and the 20% of the population with the highest 
values of this ratio. 

It is reassuring that we recover a similar dependence 
of the linear bias on halo mass when labelling haloes by 
Knax/Vioo as other authors foun d using the concentration 
parameter (| Wechsler et al.l 120061 ) . This trend is the oppo- 
site to that r ecovered when halo samples ar e split by for- 
mation time. iGao. Springel, fc White! l|2005h found no de- 
pendence of the clustering signal on halo formation time 
for massive haloes. This is puzzling since formation time 
and concentratio n are correlated, albeit with sc atter (e.g. 
iNeto et all 120071 1. ICroton. Gao. fc White! (|2007h have ar- 
gued that this suggests that an as yet unknown halo prop- 
erty is a more fundamental property in terms of determin- 
ing the clustering strength (for theoretical explanations of 



the physical basis of assembly bias see e.g. IZentnerl [20071 : 
lAriel Keselman fc Nusser 1120071 ; [Dalai et al.ll2008l ). 

The second order bias parameter, 62, displays qualita- 
tively similar dependences on mass and Knax/V^oo to 61 with 
the difference that 62 is negative around 8 c /a(M, z) ~ 1. 
The third and fourth order bias parameters are more com- 
plicated, being essentially independent of mass until peaks 
S c /a(M, z) ~ 2 — 3 are reached, where there is a dip in bias 
before a rapid increase for rarer peaks. The dependence on 
V m ax/V2oo is also different at third and fourth order. 

We compared our measurements for the bias parame- 
ters with analytic predictions. For haloes corresponding to 
rare peaks, the trend in linear bias ve rsus peak height is 
interm ediate between the predictions of iMo. Jing. fc White! 
l| 19971 ). which are based on Press & Schechter's (1974) the- 
ory for the abundance of haloes and the spherical collapse 
model, and the calculation of Sheth, Mo & Tormen (2001) 
and Scoccimarro et al. (2001), based on ellipsoidal collapse 
and an improved estimate of the halo mass function. Both 
analytic calculations predict a weaker dependence of 61 on 
peak height around 5 c /a(M, z) ~ 1 than we find in the sim- 
ulation. The comparison between the simulation measure- 
ments and the analytic predictions is similar for &2- For the 
third and fourth order bias parameters, the simulation re- 
sults are in good agreement with the analytic predictions 
for modest peaks. For rare peaks, the bias parameters mea- 
sured from the similar are again in between the two analytic 
predictions. 

Observations of clustering are already entering the 
regime in which our simulation can play an important role 
in interpreting the measurements. Existing observations of 
high redshift quasar clustering suggesting that these ob- 
jects live in haloes corresponding to ~ 5 — 6 sigma peaks 
in the matter distribution at z — 4 (White, Martini & Cohn 
2007). Future galaxy surveys, due to the volume covered 
and number of galaxies targeted, will yield measurements 
of clustering with unprecedented accuracy, to higher orders 
than the two-point function. The measurements presented 
in this paper will provide invaluable input to future models 
of galaxy clustering based on halo occupation distribution 
models, which have been modified such that galaxy cluster- 
ing is a function of mass and a second halo property. 
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